18.338J/16.394J: The Mathematics of Infinite Random Matrices 
Numerical Methods in Random Matrices 



Per-Olof Persson 
Handout #7, Tuesday, September 28, 2004 



1 Largest Eigenvalue Distributions 

In this section, the distributions of the largest eigenvalue of matrices in the /3-ensembles are studied. His- 
tograms arc created first by simulation, then by solving the Painleve II nonlinear differential equation. 

1.1 Simulation 

The Gaussian Unitary Ensemble (GUE) is defined as the Hermitian n x n matrices A, where the diagonal 
elements Xjj and the upper triangular elements Xjk = Ujk + iVjk are independent Gaussians with zero-mean, 
and 

Vax(xjj) = l, l<j<n, 
Var(ujfc) = Var(vjfc) = ^ 1 < j < k < n. 

Since a sum of Gaussians is a new Gaussian, an instance of these matrices can be created conveniently in 
MATLAB: 

A=randn(n)+i*randn(n) ; 
A=(A+A')/2; 

The largest eigenvalue of this matrix is about 2y / n. To get a distribution that converges as n — > oo, the 
shifted and scaled largest eigenvalue A( nax is calculated as 

A raax = n ^ ( A max - 2y/fl) . (2) 

It is now straight-forward to compute the distribution for by simulation: 

for ii=l: trials 

A=randn(n)+i*randn(n) ; 
A=(A+A')/2; 
lmax=max(eig(A)) ; 

lmaxscaled=n" (1/6) * (lmax-2*sqrt (n) ) ; 
'/, Store lmax 
end 

7, Create and plot histogram 

The problem with this technique is that the computational requirements and the memory requirements grow 
fast with n, which should be as large as possible in order to be a good approximation of infinity. Just 
storing the matrix A requires n 2 double-precision numbers, so on most computers today n has to be less 
than 10 . Furthermore, computing all the eigenvalues of a full Hermitian matrix requires a computing time 
proportional to n 3 . This means that it will take many days to create a smooth histogram by simulation, 
even for relatively small values of n. 



To improve upon this situation, another matrix can be studied that has the same eigenvalue distribution 
as A above. In [1], it was shown that this is true for the following symmetric matrix when (5 = 2: 



/7V(0,2) X (n-i)0 

X(n-l)P N(0,2) X(n-2)0 



\ 



X2/3 N(0, 2) xp 

X& N(0,2)J 



(3) 



Here, N(0, 2) is a zero-mean Gaussian with variance 2, and Xd is the square-root of a x 2 distributed number 
with d degrees of freedom. Note that the matrix is symmetric, so the subdiagonal and the superdiagonal are 
always equal. 

This matrix has a tridiagonal sparsity structure, and only 2n double-precision numbers are required to 
store an instance of it. The time for computing the largest eigenvalue is proportional to n, cither using 
Krylov subspace based methods or the method of bisection [2] . 

In MATLAB, the built-in function eigs can be used, although that requires dealing with the sparse 
matrix structure. There is also a large amount of overhead in this function, which results in a relatively poor 
performance. Instead, the function maxeig is used below to compute the eigenvalues. This is not a built-in 
function in MATLAB, but it can be downloaded from http://www.mit.edu/~persson/mltrid. It is based 
on the method of bisection, and requires just two ordinary MATLAB vectors as input, corresponding to the 
diagonal and the subdiagonal. 

It also turns out that only the first lOni components of the eigenvector corresponding to the largest 
eigenvalue are significantly greater than zero. Therefore, the upper-left n cu toff by ~i C utoff submatrix has the 
same largest eigenvalue (or at least very close), where 



n cu toff ~ 10ns. 



(4) 



Matrices of size 



10 12 can then easily be used since the computations can be done on a matrix of size 



only lOna = 10. Also, for these large values of n the approximation x n ~ n 1S accurate. 
A histogram of the distribution for n = 10 9 can now be created using the code below. 



n=le9; 

nrep=le4; 

beta=2; 



cutoff=round(10*n~(l/3)) ; 
dl=sqrt(n-l:-l:n+l-cutoff) '/2/sqrt(n) ; 

ls=zeros (1 ,nrep) ; 
for ii=l:nrep 

dO=randn( cutoff , 1) /sqrt (n*beta) ; 

ls(ii)=maxeig(dO,dl) ; 
end 



ls=(ls-l)*n~(2/3)*2; 



histdistr (Is ,-7:0.2:3) 



where the function histdistr below is used to histogram the data. It assumes that the histogram boxes are 
equidistant. 

function [xmid , H] =histdistr (Is , x) 



dx=x(2)-x(l) ; 
H=histc (Is ,x) ; 
H=H(l:end-l) ; 




H=H/sum(H)/dx; 

xmid=(x(l:end-l)+x(2: end) ) / 2 ; 

bar (xmid.H) 
grid on 

The resulting distribution is shown in Figure 1, together with distributions for (3 = 1 and (3 — 4. The 
plots also contain solid curves representing the "true solutions" (see next section). 

1.2 Painleve II 

Instead of using simulation to plot the distributions of the largest eigenvalues, it can be computed from the 
solution of the Painleve II nonlinear differential equation [3] : 

q" = sq + 2q 3 (5) 

with the boundary condition 

q(s) r~j Ai(s), as s -> oo. (6) 
The probability distribution /2(s), corresponding to (3 = 2, is then given by 

Ma) = ±F 2 (a), (7) 

where 

(x - s)q(x) 2 dx ) . (8) 



The distributions for (3 = 1 and (3 = 4 are the derivatives of 

Fi(s) 2 = F 2 {s)e-^^ x)dx (9) 

and 

/ \2 / 1 J-°° q ( x ) dx , -i /~ q(x)dx \ 2 

F,U) -*(.) : . HO) 



To solve this numerically using MATLAB, first rewrite (5) as a first-order system: 



d_ fq 

ds W 



q 

sq + 2q :: 



(11) 



This can be solved as an initial-value problem starting at s = Sq = sufficiently large positive number, and 
integrating along the negative s-axis. The boundary condition (6) then becomes the initial values 



q'(so) 



Ai(so) 
Ai'(so). 



(12) 



Although the distributions can be computed from q(s) as a post-processing step, it is most convenient 
to add a few variables and equations to the ODE system. When computing ^(s), the quantity I(s) = 
J (x — s)q{x) 2 dx is required. Differentiate this twice to get 



I'(s) 



q(x) 2 dx 



(13) 



and 



I"(s)= q(s) s 



(14) 



Add these equations and the variables I(s),I'(s) to the ODE system, and the solver will do the integration. 
This is not only easier and gives less code, it will also give a much more accurate solution since the same 
tolerance requirements are imposed on J(s) as on the solution q(s). 

In a similar way, the quantity J(s) = q(x) dx is needed when computing -Fi(s) and F^s). This is 
handled by adding the variable J(s) and the equation J'(s) = —q(s). 

The final system now has the form 



with the initial condition 



d_ 

ds 



q'(so) 
I(s ) 

I'M 



fq\ 

q' 
I 

r 



i q' \ 

sq + 2q 3 

r 

q 2 

V -q J 



(15) 



Ai(« ) 
Ai'(so) 
f™(x-s )M(xydx 

Ai(s Q ) 2 
V J~M(z)dx J 



(16) 



This problem can be solved in just a few lines of MATLAB code using the built-in Runge-Kutta based ODE 
solver ode45. First define the system of equations as an inline function 

deq=inline(' [y(2) ; s*y(l)+2*y(l)"3; y(4) ; y(l)~2; -y (1)] ' , ' s ' , 'y ' ) ; 

Next specify the integration interval and the desired output times. 

s0=5; 
sn=-8; 

sspan=linspace(sO , sn, 1000) ; 

The initial values can be computed as 

y0= [airy (s0) ; airy(l.sO); ... 

quadl (inline (' (x-sO) .*airy(x) . ~2 ' , 'x ' , ' sO ' ) , sO ,20 , le-25,0, sO) ; ... 
airy(s0)~2; quadl (inline (' airy (x) ») , sO ,20 , le-18)] ; 



where the quadl function is used to numerically approximate the integrals in (16). Now, the integration 
tolerances can be set and the system integrated: 

opts=odeset ( 'reltol ' , le-13 , 'abstol ' , le-15) ; 
[s ,y] =ode45 (deq, sspan,yO ,opts) ; 

The five dependent variables are now in the columns of the MATLAB variable y. Using these, F 2 (s), -Pi(s), 
and i*4(s) become 

F2=exp(-y(:,3)); 
Fl=sqrt(F2.*exp(-y(: ,5))) ; 

F4=sqrt(F2) .*(exp(y(: ,5)/2)+exp(-y( : ,5)/2))/2; 
s4=s/2"(2/3) ; 

The probability distributions f 2 {s), fi(s), and /4(s) could be computed by numerical differentiation: 

f 2=gradient(F2,s) ; 
f l=gradient(Fl,s) ; 
f 4=gradient (F4 , s4) ; 

but it is more accurate to first do the differentiation symbolically: 

f 2 (s) = -I'(s)F 2 (s) (17) 

^ = opTT (f2(s) + q(s)F 2 (s))e- J M (18) 
2F 1 (s) 

U(s) = ^ (h(s) (2 + e J W + e~ J ^) + F 2 (s)q(s) (e~ J ^ - e J <»>)) (19) 

and evaluate these expressions: 
f2=-y(: ,4) .*F2; 

f 1=1/2. /Fl.*(f2+y( : ,1) .*F2) .*exp(-y( : ,5)) ; 
f4=l/2~(l/3)/4./F4.*(f2.*(2+exp(y(: , 5) )+exp(-y ( : ,5)))+ ... 

F2.*y(: ,1) .*(exp(-y(: ,5))-exp(y(: ,5)))) ; 

Finally, plot the curves: 

plot(s,fl,s,f2,s4,f4) 
legend('\beta=l' , '\beta=2' , '\beta=4') 
xlabeK ' s ' ) 

ylabel('f_\beta(s) ' , 'rotation' ,0) 
grid on 

The result can be seen in Figure 2. 




2: The probability distributions /i(s), f2(s), and f±(s), computed using the Painlevc II solution. 



2 Eigenvalue Spacings Distributions 



Another quantity with an interesting probability distribution is the spacings of the eigenvalues of random 
matrices. It turns out that the eigenvalues are almost uniformly distributed, which means that every random 
matrix gives a large number of spacings. The distributions can then be efficiently computed by simulation. 

Two other methods are used to compute the spacings distribution - the solution of the Painleve V 
nonlinear differential equation and the eigenvalues of the Prolate matrix. Finally, the results are compared 
with the spacings of the zeros along the critical line of the Riemann zeta function. 

2.1 Simulation 

As before, the simulations are made with matrices from the Gaussian Unitary Ensemble. The normalized 
spacings of the eigenvalues Ai < A2 < . . . < Xn are computed according to 



where (3 = 2 for the GUE. The distribution of the eigenvalues is almost uniform, with a slight deviation 
at the two ends of the spectrum. Therefore, only half of the eigenvalues are used, and one quarter of the 
eigenvalues at each end is discarded. 

Again, to allow for a more efficient simulation, the tridiagonal matrix (3) is used instead of the full 
Hcrmitian matrix. In this case, all the eigenvalues arc computed, which can be done in a time proportional 
to n 2 . While this could in principle be done using the MATLAB sparse matrix structure and the eigs 
function, the more efficient trideig function is used below to compute all the eigenvalues of a symmetric 
tridiagonal matrix. It can be downloaded from http://www.mit.edu/~persson/mltrid. 

The histogram can now be computed by simulation with the following lines of code. Note that the 
function chi2rnd from the Statistics Toolbox is required. 

n=1000; 
nrep=1000; 
beta=2 ; 

ds=zeros (1 ,nrep*n/2) ; 
for ii=l:nrep 

l=trideig(randn(n,l) , sqrt (chi2rnd( (n-1 :-l : 1) '*beta)/2)) ; 
d=dif f (1 (n/4 : 3*n/4) ) /beta/pi . *sqrt (2*beta*n-l (n/4 : 3*n/4- 1 ) . ~2) ; 
ds ( (ii-1) *n/2+l : ii*n/2)=d; 
end 

histdistr (ds , : . 05 : 5) ; 

The resulting histogram can be found in Figure 3. The figure also shows the expected curve as a solid line. 
2.2 Painleve V 

The probability distribution p(s) for the eigenvalue spacings when (3 = 2 can be computed with the solution 
to the Painleve V nonlinear differential equation (see [4] for details): 




k ps n/2 



(20) 




(21) 



with the boundary condition 




(22) 



Then p(s) is given by 



p(s) = ^E(s) 



(23) 




Figure 3: Probability distribution of consecutive spacings of random matrix eigenvalues (1000 repetitions, 
n = 1000) 



where 

E(s) = exp (J ~Y' dt J ■ ( 24 ) 

Explicit differentiation gives 

p(s) = (nsa'iTTs) - a(ns) + o(-ks) 2 ) E(s). (25) 
The second-order differential equation (21) can be written as a first-order system of differential equations: 

(26) 



dt \<T'J \-ly/{<T-t<T')(t<T' -<7+(<7') 2 ) 

This is solved as an initial-value problem starting at t = to = very small positive number. The value t = 
has to be avoided because of the division by t in the system of equations. This is not a problem, since the 
boundary condition (22) provides an accurate value for a(to) (as well as E{to/n)). The boundary conditions 
for the system (26) then become 

To be able to compute E(s) using (24), the variable 

I(t) = f*^T dt ' C 28 ) 
is added to the system, as well as the equation 4:1 = j. The corresponding initial value is 

7 ~ 2^2 







-) 






7T 2 J 



Putting it all together, the final system is 









-( 







fV^-toO^'-a + ^O 2 ) (30) 



with boundary condition 




_ W I _ ^ 

7T 7T 
\ 7T 2jT 2 

This system is defined as an inline function in MATLAB: 

desig=inline( [' [y(2) ; -2/t*sqrt ( (y (l)-t*y (2) ) * ' ... 

' (t*y(2)-y(l)+y(2)~2)); y (1) /t] '] , ' t ' , 'y ' ) ; 

Specify the integration interval and the desired output times: 

t0=le-12; 
tn=16; 

tspan=linspace(tO,tn, 1000) ; 
Set the initial condition: 

y0= [-t0/pi- (tO/pi) ~2; -l/pi-2*t0/pi ; -t0/pi-t0~2/2/pi~2] ; 



2\ 



(31) 
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Figure 4: Painlcvc V (left), E(s) and p(s) (right). 

Finally, set the integration tolerances and call the solver: 

opts=odeset('reltol' ,le-13, 'abstol' , le-14) ; 
[t ,y] =ode45(desig, tspan,yO , opts) ; 

The solution components are now in the columns of y. Use these to evaluate E(s) and p(s): 

s=t/pi ; 

E=exp(y( : ,3)) ; 

p=l./s.~2.*E.*(t.*y(:,2)-y(:,l)+y(:,l)."2); 
p(l)=2*s(l); '/, Fix due to cancellation 

A plot of p(s) can be made with the command 

plot(s,p) 
grid on 

and it can be seen in Figure 4. Plots are also shown of E(s) and <r(i). 
2.3 The Prolate Matrix 

Another method to calculate the distribution of the eigenvalue spacings is to compute the eigenvalues A; of 
the operator 

f(y) - f Q(x, y)f(y) dy, Q(x, y) = sin p ~ ^ . (32) 
J-i (x-y)TT 

Then E(2t) = (1 — A^) , and p(s) can be computed as before. To do this, first define the infinite symmetric 
Prolate matrix: 



Ac 



( a-o ai 



(33) 



with ao = 2w, = (sin2irwk)/irk for k = 1, 2, . . ., and < w < ^. A discretization of Q(x, y) is achieved 
by setting w — t/n and extracting the upper-left n x n submatrix A n of ^oo. 



Below, the full matrix A n is used, and all the eigenvalues are computed in n 3 time using the eig function. 
However, A n commutes with the following symmetric tridiagonal matrix [5], and therefore has the same 
eigenvectors: 



/ Oil Pi 
Pi a 2 



T„ 



\ 



Pn-1 Oi-n-1 Pn-1 

Pn-1 OL n J 



(34) 



where 



Oik 



n+l 



k(n — k). 



k) 2 cos 2irw 



(35) 



It is then in principle possible to use the new techniques described in [6] to compute all the eigenvalues and 
eigenvectors of T n in n 2 time, and then get the eigenvalues of A n by dot products. This is not done in this 
example. 

The code for computing E(s) is shown below. This time, p(s) is evaluated by numerical differentiation 
since no information about the derivative of E(s) is available. 

s=0:0.01:5; 
n=100; 

E0=zeros(size(s)) ; 
for ii=l : length(s) 

Q=gallery ( 'prolate ' ,n,s(ii)/2/n) ; 

E0(ii)=prod(l-eig(Q)); 
end 

pO=gradient (gradient (E0, s) ,s) ; 

To improve the accuracy in E(s), Richardson extrapolation can be used. This is done as follows, where the 
values are assumed to converge as 1/n 2 : 



'/, ... Compute s and E using Painleve V in previous section 



Es=zeros( length (t) ,0) ; 
El=zeros(size(s)) ; 
for n=20*2.-(0:3) 
for ii=l : length(s) 

Q=gallery( 'prolate ',n,s(ii)/2/n); 
El(ii)=prod(l-eig(Q)); 
end 

Es=[Es,El] ; 
end 



for ii=l:3 

max(abs(Es-E( : , ones(l , size (Es , 2) ) ) ) ) 
Es=Es(: ,2:end)+diff (Es , 1 , 2)/ (2~ (ii+1) -1) ; 

end 

max(abs (Es-E) ) 



The errors maxo< s <5 \Ei(s) — E(s)\ are shown in Table 1, for n = 20,40,80, and 160. The error after all 
extrapolations is of the same order as the "exact solution" using Painleve V. 



2.4 Riemann Zeta Zeros 



It has been observed that the zeros of the Riemann zeta function along the critical line Re(z) = i behave 
similar to the eigenvalues of random matrices in the GUE. Here, the distribution of the scaled spacings of the 



N 


Error 


Error 1 


Error 2 


Error 3 


20 


0.2244 








40 


0.0561 


0.7701 






80 


0.0140 


0.0483 


0.5486 




160 


0.0035 


0.0032 


0.0323 


2.2673 




•10- 3 


•10-' 


•10-* 


•io- 11 



Table 1: Difference between Prolate solution Ei(s) and Painleve V solution E(s) after 0, 1, 2, and 3 
Richardson extrapolations. 

zeros is compared to the corresponding distribution for eigenvalues computed using the Painleve V equation 
from the previous chapters. 

Define the nth zero 7„ = n th as 

C Q + i'7n) =0, < 7l < 72 < ■ • • (36) 

Compute a normalized spacing: 

7n 

7n = : = In ■ 

av spacing near 7„ 

Zeros of the Riemann zeta function can be downloaded from [7]. Assuming that the MATLAB variable 
gamma contains the zeros, and the variable offset the offset, these two lines compute the consecutive 
spacings 7 n +i — j n and plot the histogram: 

delta=dif f (gamma) /2/pi . *log( (gamma (1 : end-l)+of f set) /2/pi) ; 
histdistr(delta, 0:0. 05:5.0) ; 

The result can be found in Figure 5, along with the Painleve V distribution. The curves are indeed in good 
agreement, although the number of samples here is a little to low to get a perfect match. 
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Figure 5: Probability distribution of consecutive spacings of Riemann zeta zeros (30,000 zeros, n w 
10 12 ,10 21 ,10 22 ) 



